    module LU_tdiag
    integer,parameter :: DB=8

    contains
    subroutine LU_factorize(sm, mt)
    implicit none
    integer, intent(in) :: sm
    real(DB),dimension(3,sm), intent(inout) :: mt
    integer::g, i, k, itemp
    real(DB) :: mult,pivot

    do k=1, sm-1
        ! LU factorize
        if (mt(2,k) .eq. 0._8) cycle
        mult = mt(1,k+1)/mt(2,k)
        mt(1,k+1)= mult
        mt(2,k+1)= mt(2,k+1) - mult * mt(3,k)
    end do

    do i =1, sm
        if (mt(2,i) .ne. 0._8) mt(2,i) = 1/ mt(2,i)
    end do
    
    end subroutine

      
    subroutine solve_LU (sm, mt, rhs, sol)
    implicit none
    integer, intent(in) :: sm
    real(DB),dimension(3,sm), intent(in) :: mt
    real(DB),dimension(sm), intent(in) :: rhs
    real(DB),dimension(sm), intent(out) :: sol
    real(DB),dimension(sm) :: y
    real(DB) :: sum
    integer :: i, j
    
    ! forward
    y(1) = rhs(1)
    do i=2, sm
        y(i) = rhs(i)-mt(1,i)*y(i-1)
    end do
    
    ! backward
    i = sm
    sol(sm) = y(sm)*mt(2,i)
    do i=sm-1, 1, -1
        sol(i) = mt(2,i)*(y(i) - mt(3,i)*sol(i+1))
    end do
    
    end subroutine

    end module

module LU_general
    integer,parameter :: DB=8

    contains
    subroutine LU_factorize(sm, mt)
    implicit none
    integer, intent(in) :: sm
    real(DB),dimension(sm,sm), intent(inout) :: mt
    integer::g, i, j, k, itemp
    real(DB) :: mult,pivot

    do k=1, sm-1
        ! LU factorize
        do i=k+1, sm
            mult = mt(k,i)/mt(k,k)
            if (mult .eq. 0._8) cycle
            mt(k,i)= mult
            do j=k+1, sm
                mt(j,i)= mt(j,i) - mult * mt(j,k)
            end do
        end do
    end do

    do i =1, sm
        mt(i,i) = 1/ mt(i,i)
    end do
    
    end subroutine

    subroutine solve_LU (sm, mt, rhs, sol)
    implicit none
    integer, intent(in) :: sm
    real(DB),dimension(sm,sm), intent(in) :: mt
    real(DB),dimension(sm), intent(in) :: rhs
    real(DB),dimension(sm), intent(out) :: sol
    real(DB),dimension(sm) :: y
    real(DB) :: sum
    integer :: i, j
    
    ! forward
    y(1) = rhs(1)
    do i=2, sm
        y(i) = rhs(i)
        do j=1, i-1
            y(i) = y(i)-mt(j,i)*y(j)    
        end do
    end do
    
    ! backward
    i = sm
    sol(sm) = y(sm)*mt(i,i)
    do i=sm-1, 1, -1
        sol(i) = y(i)
        do j=i+1, sm
            y(i) = y(i) - mt(j,i)*sol(j)
        end do
        sol(i) = mt(i,i)*y(i)
    end do
    
    end subroutine
          


end module
